/*
 * GetCharges.c
 *
 *  Created on: Apr 13, 2010
 *      Author: marchi
 */
#include "stdio.h"
#include "stdlib.h"
#include "statutil.h"
#include "typedefs.h"
#include "smalloc.h"
#include "vec.h"
#include "copyrite.h"
#include "statutil.h"
#include "tpxio.h"
#include "fftgrid.h"
#include "index.h"
#include "pbc.h"
#include "math.h"
#include "get_charges.h"


#define VERTEX 8

const rvec cube[VERTEX]={{0,0,0},{1,0,0},
		{0,1,0},{0,0,1},{1,1,0},
		{1,0,1},{0,1,1},{1,1,1}};


static real TriDiff(real g, real cube){
	return((g-cube>=0.0)? cube-g+1 : g-cube+1);
}

static void TriLinear(real gx, real gy, real gz, real q, real *qq){
	int n;
	real dd;
	real sum=0;
	for(n=0;n<VERTEX;n++){
		dd=TriDiff(gx,cube[n][XX])*TriDiff(gy,cube[n][YY])*TriDiff(gz,cube[n][ZZ]);
		qq[n]=dd*q;
		sum+=qq[n];
	}
}

int get_charges(real*** Rho_slt, real*** Rho_sol,real* Charges,int nx,int ny,int nz,int nind, atom_id *cindex, t_trxframe* fr){

	int i;
	real** xa;
	matrix oc,co;
	int n,m;


	for(n=0;n<DIM;n++){
		for(m=0;m<DIM;m++){
			co[m][n]=fr->box[n][m];
		}
	}


	calc_recipbox(co,oc);

	int natoms=fr->natoms;
	real x1,y1,z1,r1,s1,t1,gx,gy,gz;
	int mx,my,mz;
	real qq[8];
	int ox,oy,oz,tx,ty,tz;

	xa = malloc(sizeof(real *) * natoms);
	for (i = 0; i < natoms; i++) {
		xa[i] = malloc(sizeof(real) * DIM);
	}
	for (i = 0; i < natoms; i++) {
		xa[i][XX] = oc[XX][XX] * fr->x[i][XX] + oc[XX][YY] * fr->x[i][YY] + oc[XX][ZZ] * fr->x[i][ZZ];
		xa[i][YY] = oc[YY][XX] * fr->x[i][XX] + oc[YY][YY] * fr->x[i][YY] + oc[YY][ZZ] * fr->x[i][ZZ];
		xa[i][ZZ] = oc[ZZ][XX] * fr->x[i][XX] + oc[ZZ][YY] * fr->x[i][YY] + oc[ZZ][ZZ] * fr->x[i][ZZ];
	}

	for (i = 0; i < natoms; i++) {
		x1 = xa[i][XX];
		y1 = xa[i][YY];
		z1 = xa[i][ZZ];
		r1 = (real) nx * (x1 - rint(x1 - 0.5));
		s1 = (real) ny * (y1 - rint(y1 - 0.5));
		t1 = (real) nz * (z1 - rint(z1 - 0.5));
		mx = (int) r1;
		my = (int) s1;
		mz = (int) t1;

		gx = r1 - (real) mx;
		gy = s1 - (real) my;
		gz = t1 - (real) mz;
		TriLinear(gx, gy, gz, Charges[i], qq);
		for (n = 0; n < VERTEX; n++) {
			tx = (int) cube[n][XX];
			ty = (int) cube[n][YY];
			tz = (int) cube[n][ZZ];
			ox = (nx > mx + tx) ? mx + tx : mx + tx - nx;
			oy = (ny > my + ty) ? my + ty : my + ty - ny;
			oz = (nz > mz + tz) ? mz + tz : mz + tz - nz;
			Rho_sol[ox][oy][oz] += qq[n];
		}
	}
	for (int ii=0; ii < nind;ii++) {
		i=cindex[ii];
		x1 = xa[i][XX];
		y1 = xa[i][YY];
		z1 = xa[i][ZZ];
		r1 = (real) nx * (x1 - rint(x1 - 0.5));
		s1 = (real) ny * (y1 - rint(y1 - 0.5));
		t1 = (real) nz * (z1 - rint(z1 - 0.5));
		mx = (int) r1;
		my = (int) s1;
		mz = (int) t1;
		gx = r1 - (real) mx;
		gy = s1 - (real) my;
		gz = t1 - (real) mz;
		TriLinear(gx, gy, gz, Charges[i], qq);
		for (n = 0; n < VERTEX; n++) {
			tx = (int) cube[n][XX];
			ty = (int) cube[n][YY];
			tz = (int) cube[n][ZZ];
			ox = (nx > mx + tx) ? mx + tx : mx + tx - nx;
			oy = (ny > my + ty) ? my + ty : my + ty - ny;
			oz = (nz > mz + tz) ? mz + tz : mz + tz - nz;
			Rho_sol[ox][oy][oz] -= qq[n];
			Rho_slt[ox][oy][oz] += qq[n];
		}
 	}
	for (i = 0; i < natoms; i++) {
		free(xa[i]);
	}
	free(xa);

	return 0;
}
